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CNl ■ The hydrodynamic description of probabilistic ballistic annihilation, for which no conservation 

laws hold, is an intricate problem with hard sphere-like dynamics for which no exact solution exists. 
We consequently focus on simplified approaches, the Maxwell and very hard particles (VHP) models, 
which allows us to compute analytically upper and lower bounds for several quantities. The purpose 
is to test the possibility of describing such a far from equilibrium dynamics with simplified kinetic 
models. The motivation is also in turn to assess the relevance of some singular features appearing 
within the original model and the approximations invoked to study it. The scaling exponents are 
first obtained from the (simplified) Boltzmann equation, and are confronted against Monte Carlo 
simulation (DSMC technique). Then, the Chapman-Enskog method is used to obtain constitutive 
relations and transport coefficients. The corresponding Navier-Stokes equations for the hydrody- 

Snamic fields are derived for both Maxwell and VHP models. We finally perform a linear stability 
analysis around the homogeneous solution, which illustrates the importance of dissipation in the 
possible development of spatial inhomogeneities. 

a ~ 

5j ' PACS numbers: 05.20.Dd,47.20.-k,45.05.+x,82.20.Nk 



INTRODUCTION 



The possibility to describe in terms of hydrodynamic equations the evolution of a system where some physical 
q i quantities are not conserved is a challenging problem of non-equilibrium statistical mechanics. Several questions have 
to be faced as for example the validity of the underlying (and in practice approximate) kinetic theory, the choice of the 
l— ~~ ', hydrodynamical fields that are supposed to describe the relevant excitations in the problem, or the consistency of the 
method itself that is used to deduce the coarse-grained description from the kinetic theory. Much attention has been 
recently paid to such questions, mainly in the field of granular gas dynamics (see e.g PIS SB)- In such systems, 
• the kinetic energy is not conserved, while the linear momentum and number of particles are. However, even for low 
dissipation, the derivation of the hydrodynamic relations, based on a hard-sphere-like Boltzmann equation is not a 
simple task and several approximations have to be invoked [2j . These difficulties lead to consider some simpler models 
by choosing ad- hoc collision term in the Boltzmann equation. The so-called Maxwell and very hard particles (VHP) 
models [5l |ol are particularly interesting and reproduce some qualitative features of the granular gas of inelastic hard 
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Another class of problems for which not only energy but also the density and momentum are not conserved is 
probabilistic ballistic annihilation (PBA). In such a system, the particles move ballistically between collisions. When 
two particles meet, they undergo an instantaneous collision and are removed from the system with probability p or 
I ' undergo an elastic scattering with probability (1 — p). Since collisions are assumed to be instantaneous, two body 
events only are taken into account. The PBA model was first introduced in one dimension in In the limit 

p — > 0, where density, momentum and kinetic energy are conserved, one recovers a system of hard spheres for which 
the hydrodynamic equations are well know n |l5l 1161 Il7l . The other limit p — 1 (pure annihilation) has been the 
object of some work p~8l Il9l I20I l2ll 123 . 123L 12 ll l'25l l26j. It was shown that in the long time limit the annihilation 
dynamics is exactly described by the Boltzmann equation in dimensions higher than one |24| . This may qualitatively 
be understood by the fact that the density of the gas decays and, at late times, the packing fraction is very low. This 
fact lead to conjecture that the Boltzmann equation is an adequate description of PBA at late times for p > |27j . 

Given that p may be considered as a perturbation parameter allowing to recover the elastic limit, the PBA model 
is particularly interesting in view of testing the relevance and validity of the hydrodynamic description in general, 
which is a controversial issue and a long term goal of the present work . The analytical treatment with usual hard 
sphere dynamics however appears to be quite involved |28| . and we study here the simplified Maxwell and VHP 
versions of PBA. The motivation is here is not only to test the ability of simplified kinetic models to mimic the hard 
sphere dynamics for a model far from equilibrium (and with no conserved quantity, a more severe situation than 
that of granular gases) but also to shed some light on some peculiar features obtained in the hydrodynamic study of 
Ref. [2£|. In particular, this work exhibited divergent transport coefficients for a critical value of p. We will see that 
such singularities are absent in the simplified approaches, which may indicate that they are not associated with any 
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physically relevant phenomenon. It will also appear that Maxwell and VHP approaches provide useful bounds for the 
hard sphere dynamics, so that similar inequalities as those found in |23t 29j concerning the scaling exponents can be 
obtained. 

The paper is organized as follows. In section ^ we introduce the Boltzmann equation for both Maxwell and 
VHP models of probabilistic ballistic annihilation, as well as the balance equations for the coarse-grained fields. In 
section IIIII we briefly describe the Chapman-Enskog scheme while section IIVI is devoted to the Maxwell model. We 
first find the homogeneous state, and solve the corresponding homogeneous balance equations. To first order in the 
Chapman-Enskog expansion we then study the effect of a small spatial inhomogeneity. We follow the traditional route 
to compute the transport coefficients, which consists in truncating the first-order velocity distribution function to the 
first nonzero term in a Sonine polynomial expansion |28|. We then show that this truncation does not constitute an 
approximation for the transport coefficients since they can be obtained by solving the Maxwell model exactly to first 
order. The VHP model is subsequently investigated in section We first find the homogeneous cooling state, and 
then solve the corresponding homogeneous equations. We implement Monte Carlo simulations in order to check the 
decay exponents found analytically. Next, we establish the transport coefficients to first order in the Chapman-Enskog 
expansion before presenting a comparison of the transport coefficients of the different models. In Sec. IVIll we finally 
perform a linear stability analysis of the Navier-Stokes hydrodynamic equations around the spatially homogeneous 
state, and compare the results with PBA of hard spheres. We show as well that the second order decay rates may 
accurately be neglected. Our main findings and conclusions are summarized in section fVIIII 

Since the underlying calculations of this paper are cumbersome, we present only the main steps in order to focus 
onto the more relevant results. Further technical details or explanations may be found in |28| and for convenience, 
Appendix lAl contains a summary of the notations used. 



II. THE BALANCE EQUATIONS 



The Boltzmann equation for the one particle distribution /(r,v;£) of particles annihilating upon collision with 
probability p reads 

(d t + vi ■ V)/(r, vi; t) = P J a [/,/] + (1 - P )J C [/,/], (1) 
where J a is the annihilation operator defined by 

J a [f,g] = -a d - 1 cp(x)v 1 T ~ x g(T,v l -t) f dv 2 vf 2 f(r,v 2 ;t) (2) 

JTR d 

and J c is the collision operator: 

Jc[f,9] = * d - 1 ^r L X I dv * v ™ [d&(b- 1 -i)g(T,Vi;t)f(w,t). (3) 

In these equations, d denotes the spatial dimension, vi 2 = |vi — V2I is the modulus of the relative velocity, Sd — 
2ir d / 2 /T(d/2) is the solid angle surface, T the Euler gamma function, vt — \/2//3m the time-dependent thermal 
velocity, (3 = (fc^T)" 1 , a is the diameter of the particles, cr is a unit vector joining the centers of two particles and 
the corresponding integral is running over the solid angle. Finally, 6 _1 an operator acting on the velocities as follows: 



b 1 v 12 = V12 - 2(vi2 • <x)<x, (4a) 
6 _1 vi = vi — (vi2 • <x)<x. (4b) 

The choice x = (x = 2) corresponds to the Maxwell (VHP) model, respectively. For hard sphere dynamics, that 
would correspond to 1 = 1, the relative velocity vi 2 gives the rate of collision and its presence makes analytical 
progress difficult. A convenient simplification to overcome this difficulty is to replace it by vf 2 v^T x where vt is 
introduced for dimensional reasons. The quantity 4>{x) which sets the relevant time scale in the problem can be freely 
chosen, and will be used in the following analysis to obtain the desired limiting behaviour in the limit p — > (see 
also for related considerations). We also note that particles interacting with an inverse power-law potential are 
described by a kinetic equation of the same form as Eq. J3J 6] . 
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In order to write hydrodynamic equations, we need to define local hydrodynamic fields: 

n(r,t) = / dvf(r,v;t), (5a) 

u(r,i) = — 1—r f dvv/(r,v;t), (5b) 
n{r, t) J Rd 

T{v ^ = < lk A i ^V 2 /(r,v;t), (5c) 
n(r, i)fc B d J R d 

where n(r,i), u(r,i), and T(r,t) are the local number density, velocity, and temperature, respectively (the latter 
definition being kinetic with no thermodynamic basis). The definition of the temperature follows from the principle 
of equipartition of energy. In Eq. (|5cl) . is the Boltzmann constant and V = v — u(r, t) is the deviation from the 
mean flow velocity. The balance equations follow from integrating the moments 1, wv, and mv 2 /2 with weight given 
by the Boltzmann equation Q . Following the same route as in |28| we thus obtain 

d t n + Vi(nttj) = -pw[f, /], (6a) 

d t Ui + —V j P ij +u j V j u i = -p-u[f,Vif], i = l,...,d, (6b) 
mn n 

2 T m 

d t T + UjVjT + ——(P^iUj + Vjqj) = p-u[f, /] - p—— -w[f, V 2 f], (6c) 

with implicit summation over repeated indices, u = (u\, . . . , ua), and 

<o\f,9] = - I d^J a [f ig }. (7) 

JB. d 

In the balance equations the pressure tensor Pjj and heat-flux q.i are defined by 

P l3 (r,t) = m [ dvViVj/Cr.vjt) = / dv/(r,v;t)Aj(V) + (8) 



ft(r,t)= / dvSi(V)/(r,v;t), (9) 



where (3 = l/(ksT) and 



DulV) = m[VlV j ~S ij ), (10) 



W) = (^-^t)k-. (ID 

As expected, when the annihilation probability p — > 0, all three coarse-grained fields n, u, and T are conserved. 

III. THE CHAPMAN-ENSKOG SOLUTION 

The Chapman-Enskog method allows from Eqs. 10 to build a closed set of equations for the hydrodynamic fields (see 
e.g 0jI3)- For this purpose, it is required to express the functional dependence of the pressure tensor and of the 
heat flux qi in terms of the hydrodynamic fields. The Chapman-Enskog approach relies on two important assumptions. 
The first one is the existence of a normal solution in which all temporal and spatial dependence of the distribution 
function /(r, v; t) may be expressed in terms of the hydrodynamic fields, /(r, v; t) = fly, n(r, t), u(r, t), T(r, t)]. The 
discussion of the relevance of this first assumption can be found elsewhere (e.g in 0). The second assumption is 
based on the separation of the microscopic time scale (the average collision time between particles) and macroscopic 
time scale (the evolution of the hydrodynamic fields and their inhomogeneities). This separation implies that the 
hydrodynamic fields are only weakly inhomogeneous, which allows for a series expansion in the gradients of the fields, 
/ = / (0) + + e 2 / (2) + . . ., where each power of the formal small parameter s is associated to a given order in 
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spatial gradients. The Chapman-Enskog method assumes the existence of an associated time derivative hierarchy: 
d/dt = d^°'/dt + ed^ /dt + e 2 d^ jdt + . . .. The insertion of these expansions in the Boltzmann equation yields 



E £ % +Vl - v )5>' /(0=pJa 



,k>0 



l>0 



E^,E^ 



l>0 



l>0 



+ (1-P)J C 



(0 



l>0 



l>0 



(12) 



The Chapman-Enskog solution is obtained upon solving the equations order by order in e. 



IV. THE MAXWELL MODEL 
A. The homogeneous state 

To zeroth order in the gradients, Eq. I|12H gives 

^ (0) / (0) =pJa[/ (0) ,/ (0) ] + (1 -P)Jc[/ (0) ,/ (0) ]. (13) 

This equation has a solution, describing the homogeneous state, and which obeys the scaling relation 

/(°)(r,v;;) = ^/(c), (14) 

where Vt = [2/ {(3m)] 1 ' 2 is the time dependent thermal velocity, and c = V/vt, V = v — u. The existence of a scaling 
solution of the form (|14|) seems to be a general feature present in different but related contexts jToL l24l |29| . This 
solution being isotropic, one has u = 0. 

Santos and Brey 30] showed that there exists a relationship between the homogeneous solutions of the Maxwell 
model with p = and p ^ 0. We shall here briefly reproduce their arguments. It is possible to rewrite the Boltzmann 
equation 113|) for x = under the form 

4 0) / (0) (v;i') = -(C5 + C fl )n(t')/ (0) (v;t')+ / dvi /^ x (£)/ (0 W)/ (0) (vi;O, (15) 

where t' = (l—p)t, Cs = J derx(er), x(<x) = a d ~ 1 <p(x = 1)vt/S<i, and Cr = J dax{&)p/ (1 — P) is the removal collision 
frequency. Integrating Eq. (|15|l over v, the evolution of the number density is governed by dt>n(t') = —Cnn 2 (t'), 

the solution being n(t') — no/(l + n^Cat'), where uq = n(t' = 0). If we define r(t') = L dsn{s)/fiQ and F(v;r) = 
/(°)(v; t')no/n(t'), then F(y; r) satisfies the Boltzmann equation without annihilation (i.e., Cr — 0). F(v; r) therefore 
evolves towards a Maxwellian, and so does / < -°- ) : we have /(c) = e~ c j-K d l 2 . 



B. The zeroth-order Chapman-Enskog solution 

Since is isotropic, to zeroth order the pressure tensor JSJ) becomes pjj — p^Sij, where p^ = nksT is the 
hydrostatic pressure, and the heat flux © becomes q^ - 1 = 0. The balance equations to zeroth order read 

d[ 0) n = -pn^\ (16a) 
a t (0) Wl = -pert®, (16b) 

(16c) 



where the decay rates are 











t(0) _ 
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,/ (0) 



(17a) 
(17b) 
(17c) 
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For antisymmetry reasons, one sees from Eq. <|17b(l that = 0. The calculation of £^ and £^ are straightforward 

and give = na d ~ 1 (j) M VT and £j? = 0. We have written <f> M for <p(x = 0). The temperature of the Maxwell model 
is therefore conserved in the scaling regime (time independent thermal velocity vt)- In addition, one has 

n H (t) = ^ — , (18) 

where the subscript H denotes a quantity evaluated in the homogeneous state, and $i\o) is the decay rate for t = 0. 
Note that Eq. (|18|1 was already established in Sec. IIV Al 

C. The first-order Chapman-Enskog solution 

To first order in the gradients, the Boltzmann equation ljl2|l reads 

[9 t (0) + j]/w^-[a t (1) + v 1 .v]/(°), (19) 

the operator J being defined by Eqs. (|A6(I and (|A7|) . The balance equations © to first order become 

d[ 1] n + Vi(nu<) = -pn&\ (20a) 

d\ 1) u i + —V i {nT)+u j V j u i = -wiT&, i = l,...,d, (20b) 
ran 

d[ 1] T + urf.T + jTViUi = - P T&\ (20c) 

where the decay rates are given by 

£« = -u,[fW,f(% (21a) 
n 

e« = -l,[/(»),i/,/W] + -,[/( 1 ),y 1 / (0) ], i = i ) ... > d > (2ib) 

roi>r n«T 

= _!„[/(«», /(!)] + ™ w[/ (0) )V2/ (i) ]+ - w [/(i),V/P)]. (21c) 
7i ukbI a ukb! a 

By definition is of first order in the gradients of the hydrodynamic fields; for a low density gas [j| 

/(!) = InT + B i V l In n + C, ,V ,//,. (22) 
The coefficients Ai, Bi, and depend on the fields n, V, and T. 

1. The approximate first-order Chapman-Enskog solution 

The hydrodynamic description of the flow requires the knowledge of transport coefficients, which may be determined 
from a Soninc polynomial expansion of the first order distribution function. In addition, the pressure tensor may be 
put in the form 



Pij(r, t) = p (0) % - 77 ^ViUj + VjUi - ^SijVkUk) - QSi-jVkUk, 



(23) 



where p^ = nkgT is the ideal gas pressure, r\ is the shear viscosity, and £ is the bulk viscosity which vanishes for a 
low density gas [2|. Fourier's linear law for heat conduction is 

q t = -nV,T~nV in , (24) 



where k is the thermal conductivity and /1 a transport coefficient that has no analogue in the elastic case |3lL 13/ 
The identification of Eq. (|23|l with Eq. (jHJ using the result of the first order calculation yields 



PW= [ dvAi(V)/W. (25) 
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Similarly, the identification of Eq. (|24() with Eq. © using the first order calculation leads to 

JR d 



The calculation follows the same route as in 281, and we obtain 



- JL 
m 

K 
K 

nfj, 
Tk q 



d- 1 1 



0. 



(26) 

(27a) 

(27b) 
(27c) 



where the thermal conductivity kq and shear viscosity 770 coefficients for hard spheres (used here to obtain dimen 
sionless quantities) are given by Eqs. (|A2|) and (|A3|) . respectively j^" 
given by 



The dimensionless coefficients v* and u* are 



1 j^ d dvSi{v)JAi 1 J Rd dv Sj(y)nAi 

"V J Rd dVS l (V)A 1 ' 
f J Rd dVA,(V)fiC, 



vo J Rd dVS i (V)A l 

1 / Rd dVA 3 -(V)JCi ; 

u / R „dV £><,•( V)Cy % J Rd dVD tJ (V)C 



(28a) 
(28b) 



where vq = p^/rjo, with = nksT. Note that the above relations are still exact within the Chapman- Enskog 
expansion. The approximation consists in truncating the function /W to the first nonzero term in a Sonine polynomial 
expansion: 



A(V) 
B(V) 
C(V) 



oiM(V)S(V), 
6iMCV)S(V), 
c A4(V)D(V), 



(29a) 
(29b) 
(29c) 



where ai, 61, and Co are the coefficients of the development, and A4(V) = n/ (v^n d / 2 ) exp(—V 2 /v^) is the Maxwellian 
in the homogeneous state. This allows one to compute the relations l|28|l . and one finds (see Appendix 



V2T(d/2) 

' J 4^(d-l)/2 

M V2r(d/2) 

4 7r (rf-i)/2 



rf + 2 



+ (l-p) 



d + 2 



(30a) 
(30b) 



The parameter </> M governing the collision frequency may be freely chosen to allow for a relevant comparison with 
hard sphere dynamics (see e.g. 0). We choose </> such that the transport coefficients are normalized to one for p — > 0, 
that is when all collisions are elastic. It is remarkable that for the Maxwell model a single parameter such as (f> is 
sufficient to ensure normalization of all the transport coefficients (this will not be the case in the VHP approach). 
This leads to 



iM 



47r(<*-i)/2 
V2T{d/2) ' 



(31) 



The above value turns out to be the same as the one obtained from the elastic limit of the Maxwell model of 
granular gases [T(il |. In the latter case, <fi was chosen matching the temperature decay rate with that characterizing 
the homogeneous cooling state of inelastic hard spheres. With the choice (|31|l the transport coefficients (|27|l become 

(32a) 

(32b) 
(32c) 



p 4±2 + (l- P y 
1 

0, 
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Following the same route as in [28| , the first-order distribution function (|22[) reads 

(33) 



/M(r,V;i) = -^7W(V) 
n 



Z77? n 

jSiCVJ/sViT+lDyCVJVjtii 



,2. TTie esacf first-order Chapman- Enskog solution 

By construction of the Chapman-Enskog method, the velocity moments of / are given by those of the local equi- 
librium distribution f(°>. It is then easy to show that the decay rates to first order l|21[l are equal to zero [therefore 
fi/W = 0, where the operator f2 is defined in Appendix ID] . Proceeding in a similar way as in we obtain in 
Appendix [D] the exact transport coefficients for the Maxwell model, i.e. without any approximation on the form of 
/W. This may be done by integrating the Boltzmann equation I|D1|) over V with weight mViVj and mV 2 Vi/2. With 
the choice for <f> M given by Eq. (|31|) . one finds the same transport coefficients as those given by Eqs. (|32[1 . This means 
that the truncation of f^- 1 ' to its first nonzero term in a Sonine polynomial expansion is a harmless approximation 
when looking at the transport coefficients (this is a peculiarity of the Maxwell model). In fact, it turns out that the 
transport coefficients depend only on the first term in the Sonine polynomial expansion of [T7j. For example, the 
heat current 10 may be rewritten under the form 0, 

tf^-^Z^^iT + hVtn), (34) 
2 mp 6 

where the first nonzero coefficients (a%, bi) (that may depend on n and T) in the Sonine expansion are defined by 
Eqs. (|29|l . Therefore the latter coefficients always give an exact result for the transport coefficients, but the problem 
at hand is to calculate them exactly. This turns out to be possible within the Maxwell model. 



D. Hydrodynamic equations 

Since the pressure tensor and the heat flux defined by Eqs. I|23|l and l|24|> . respectively, are of order 1 in the gradients, 
their insertion in the balance equations 10 yields contributions of order 2. Knowledge of the second order velocity 
distribution is therefore required in order to find the correct decay rates that contribute to Navier-Stokes order. 
However, we show in Sec. IVIIBl that for the Maxwell gas those contributions are equal to zero for any annihilation 
probability p. The corresponding hydrodynamic Navier-Stokes equations are given by 

dtn + V, (nm ) = -pn + #>] , (35a) 

d t Ui + — V j P ij +u j V j u i = -pv T [Z$ + i = l,...,d, (35b) 

d t T + mViT + — — {P tJ V t u 3 + Vm) = - P T[^ ] + . (35c) 

Pij and qj are given by Eqs. (|23|l with £ — 0, and ( | 24| respectively. The rates £n \ , and ma y be calculated 
using their definition l(2U|l and the distribution 28]. We find that all decay rates are equal to zero, except 

e = ^H>. (36) 
We thus have a closed set of equations for the hydrodynamic fields to the Navier-Stokes order. 



V. THE VHP MODEL 



A. The homogeneous cooling state 

Integrating the Boltzmann equation over V for x = 2, one obtains 

dn 



tU -pw(t)n, (37) 
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where 



u(t)=n(t)v T (t)a d - 1 ^ til '(c{ 2 }, 



(38) 



and (g(ci,c 2 )) = J^d dc\dc% g(cx, c 2 )/(ci)/(c 2 ) denotes the average of a function <7(ci,c 2 ) in the homogeneous 
cooling state (HCS). We have written cj) VHP for (j)(x = 2). Following the same route as in |3 113 the Boltzmann 
equation may be rewritten in the form 



(C? 2 > 



1 — a„ 



ci 



fU-0 - /«■,) / rfc 2 c? 2 /(c 2 )-^^ 



where 



and 



J R 2d dc!rfc 2 / cf 2 cf/(ci)/(c 2 ) 



/ R d dcc 2 f(c) J u2d dcidc 2 /derc 2 2 /(ci)/(c 2 ) 
![/,/] = / dd / d^c^b" 1 - l)/( Cl )/(c 2 ). 



(c?)(cf 2 



(39) 



(40) 



(41) 



The limi t ci — > of the Boltzmann equation (|39|l encodes a useful information for ballistically controlled dynam- 
ics [Hilllalll: 



(ci 2 ) [i + d 



1 — a,. 



/(0) - /(0)(c 2 ) 



1 -p 1 



(42) 



2 / p Sjci^o 

Next, we consider the first nonzero correction to the Maxwellian in a Sonine polynomial expansion of the HCS: 

f(c)=M(c)[l + a 2 S 2 (c 2 )], (43) 

where M(c) = 7r _<i / 2 e _c2 is the Maxwellian and S 2 (c 2 ) = c 4 /2 - (d + 2)c 2 /2 + d(d + 2)/8 the second Sonine poly- 
nomial [T3I • Eqs. I|42|l and (|40|l form a system of two equations for the two unknown a e and a 2 . Making use of the 
relations l|C3|l , it is a straightforward task to compute the limit in the right hand side of Eq. I|42fl [35| , which gives 



~ ~ ~ S d d 2 (d + 2) 

hm /[/,/] = -a 2 -^^— . 



Using Eq. 1)43 1). one easily obtains from Eq. l)4*0j l 



d+1 d+ 2 
a 2 - 



(44) 



(45) 



d 2d 

Note that Eqs. (|44|) and (|45|l are exact relations for which all nonlinear contributions in a 2 were kept. However, those 
nonlinear terms cancel out in each case. Making use of (c 2 2 ) = d, the insertion of Eqs. 1)44)1 and 1)43)1 in (|42|l gives 



1 - a 2 - 



d + 2 



1 + a 2 



d(d + 2) 



1 + a 2 



d(d + 2) 1 



(46) 



Eq. (|46|l admits two solutions, the first one being a 2 = and the second one a 2 = — 2[d + p(4 — d)]/[d(d + 2)p]. The 
second solution is not physical since it diverges for p = 0. Therefore a 2 = and the HCS of the VHP model within 
the approximation 1)43)1 is described by the local Maxwellian M.(c) = ir~ d / 2 e~ c . We also note that upon discussing 
the potential ambiguities resulting from such a linearization scheme in et 2 (as done in |35j,|36j), the same conclusion 
is reached. 



B. The zeroth-order Chapman-Enskog solution 

Proceeding in a similar way as already described, we obtain a set of equations formally identical to Eqs. 1)16)1 

and l|l 7)1. The calculation of the decay rates gives ^ 
therefore given by 



0, ^ 



r d-UVHP 



v T , and £i 0) = Cr'd- The HCS is 

(47a) 
(47b) 



(0). 



n H (t) = n Q (l+pt/t ) 7 ", 
T H (t)=T (l+pt/t )-^, 
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where the decay exponents are 7 n = (0)<o, 7t = £t (O)ioi an d the relaxation time i 1 = £n°^(0)+£y /2. ln other 
words, we have 

2d , 2 

7„ = — , and 7t = — : • (48) 

' 2d+V ' 2d+l v ' 

These quantities do not depend either on <fr nor on the annihilation probability p. The former result is an exact 
property of the dynamics under study (the factor </> may be absorbed into a rescaling of time t, leaving scaling 
exponents unaffected) while the latter may a priori be an artifact of the approximations made (it will however be 
shown below that the p dependence -if any- is extremely weak). If we define the root-mean-square velocity by 
v = yl (v 2 ), then from the definition (|5cl) of the temperature v(t) oc T^/ 2 (t), and from Eq. (|47bll we have v ~ t~ lv 
for long times, with j v — 7t/2. The decay exponents 7j>, and 7„, as well as the decay exponents for the Maxwell 
model, agree with the prediction of Krapivsky and Sire |23| , and satisfy the scaling constraint 7„ + 7„ = 1 , which 
essentially expresses the unicity of the relevant time scale in the problem. Moreover, making use of the expression 
for the decay ex pon ents of PBA of hard spheres 7^ and 7^ obtained to linear order in 02 and which are recalled in 
Appendix [5] p^feSj , it is easy to verify explicitly that the Maxwell and VHP models provide bounds [22j 

^T1<7^)<1, 0<7„ S (,)<^ T , (49) 



for all p £ [0, 1]. We emphasized however that the previous inequality have the status of "empirical" observations, 
and could not be anticipated from rigorous arguments. 

We performed Direct Monte Carlo Simulations (DSMC) in order to verify the decay exponents of the VHP model. 
The algorithm is similar to the one described in |27l |29|. For the sake of completeness, we briefly outline the main 
steps of the algorithm. We choose at random two different particles {i, j}. The time is then increased by VT/(N 2 vfj) 
where TV is the number of remaining particles. With probability p the two particles are removed from the system, 
and with probability 1—p their velocities are modified according to Eqs. (0J). As the fluctuations increase for small 
N , it is necessary to average over several independent realizations in order to diminish the noise. A log-log plot of the 
density n/no and the root-mean-squared velocity v/vq as a function of time gives the decay exponents (see Fig.lT)l. 
The DSMC results are in excellent agreement with the analytical predictions and the expected power-law behaviors 
are observed over several decades (see. Fig. |2J). 



C. The approximate first-order Chapman-Enskog solution 

The procedure is similar to the one followed within the Maxwell model of Sec. ITVCl (or 0), and we find 

1 



V 



K 



1^(0)*' 



d - 1 2K - 2pCn 



(0), 



X 



where X = v*[2»; - 2p&> - Sp^l +K^*{-4^ + 3p[£T + 2$"*]}, 



M* = 2 P 



d-H { v* 



(0)* 



;(0)*i 



vo J Rd dVS t (V)B t % J Rd dVS t (V)Bi ' 



(50a) 

(50b) 
(50c) 

(51) 



and f„ = £n /vq, £t — £t l v o- Truncating the function to the first term in a Sonine polynomial expansion 
as it was the case for Eqs. H29JI . the coefficients v*, v*, and j/* may be calculated with the help of Appendix |0 We 
find 



/VHP 



1 VHP 



V2T(d/2) 
47r (d-i)/2 

V2T(d/2) 

47r (d-i)/2 



' (d + 2) 2 , h , (d + 2)(d + 4) 

P n + i 1 -P) a 



p (d+2)(d+3) + (1 p) (d-l)(d + 4) 



(52a) 
(52b) 
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FIG. 1: The decay exponents 7„ and 7„ (inset) in two dimensions for the VHP model (x = 2). The analytical predictions 
7„ = 2d/(2d + 1) = 0.8 and -y v = l/(2d + 1) = 0.2 are shown by the continuous lines while the symbols stand for the DSMC 
results (obtained from approximately 300 independent runs and 10 7 initial particles) . From the above data, it appears that the 
scaling relation j n + j v — 1 is well obeyed (the deviation from 1 does not exceed 4 x 10~ 4 ) and that the scaling exponents do 
not depend on p. Note the small y scale. 




-9 -8 -7 -6 -5 -4 -3 -2-10 1 

log 10 t 

FIG. 2: Time dependence of n and v (inset) for d = 2 and p — 0.5 on a log-log scale. The initial velocity distribution is 
Gaussian. No (resp. N) is the initial (resp. remaining) number of particles, vo = v(0) is the root-mean-square velocity at 
t = 0, whereas we write v for v(t > 0). The dashed straight line is a linear interpolation giving the decay exponent of the 
power-law, and the deviations to this law for large times is due to the low number of remaining particles. 
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The free parameter ^ VHP setting the frequency collision has a priori no reason for being the same as for the Maxwell 
model. We choose this quantity such that rf{p = 0) = 1, which means that the shear viscosity for the VHP gas is 
set for vanishing p to coincide with the shear viscosity i]o of hard spheres. This allows for a better comparison of 
the transport coefficients for the Maxwell, hard sphere, and VHP models. Other choices for </> VHP are possible. The 
condition rf (0) = 1 leads to 



,VHP 



(d + 2)(d + 4)' 



(53) 



so that 



(0)* _ 



2d 



d + A 1 
AO)* _ 1 
4t " d + ^ 



(54a) 
(54b) 



The first order distribution function reads 



3 3 



2 Til 77 

— S,(V) (nViT + fiVin) + jDij^VjUi 



(55) 



where the transport coefficients are given by Eqs (|5U|) . 



D. Hydrodynamic equations 

The decay rates to first order may be calculated using the definitions (|20|) and the distribution i|55[l 2<§|, which 
gives 

e£° = 0, (56a) 
= -vr (**^iT + M*iv«n) C (56b) 
4 1} = 0, (56c) 

where 

_ d 2 (d + 2) 2 VHP V2T(d/2) 
U ~ 8(d-l) 47r( d -D/ 2 - (&7j 

The Navier-Stokes hydrodynamic equations are thus given by Eqs. Ij35(l with the decay rates (|54J) and Ij56(l . For 
consistency, second order contributions in the gradients are again needed when evaluating the decay rates. It will be 
shown in section IVII Bl that they are small corrections to the previous relations. 



VI. COMPARISON OF THE TRANSPORT COEFFICIENTS 



We compare the transport coefficients for the Maxwell, VHP, and hard sphere models (the coefficients for the latter 
model being given in |28|]). Figs. [31 S and show 77*, n*, and as a function of the annihilation probability. Note 
that once we have chosen <f>(x = 2) such that 77* — > 1 for p — > there is no reason to expect k* — ► 1 in the same limit. 
Other choices would have been possible such as enforcing k* — > 1 when p — > 0. 

From Figures El and it first appears that Mawxell and VHP models capture the essential p dependence of 
the "hard sphere" transport coefficients. In addition, they provide in most cases lower and upper bounds for 77*, 
k* and //*. However, as already pointed out in |2q. for strong annihilation probability p ~ p^, the hard sphere 
thermal conductivity and "Fourier" coefficient fj, diverge (see Figs 01 and which leads to a violation of the VHP 
upper bound for k and fi in the vicinity of pd- The fact that VHP and Maxwell models lead to smooth and regular 
transport coefficients for all values of p gives a hint that the hard sphere divergence obtained in previous work [28| is 
a possible artifact of the underlying approximations and probably does not point towards a change of behavior nor 
a qualitative difference in the scaling or transport properties. This point will be further discussed in the concluding 
section. We finally note that an o priori similar deficiency was already reported for the Maxwell model of inelastic 
hard spheres |10| . 
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FIG. 3: Dimensionless shear viscosity rf as a function of the annihilation probability p for the Maxwell (thin continuous line), 
VHP (dashed line), and hard spheres models (thick continuous line). 
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FIG. 4: Reduced thermal conductivity k* as a function of the annihilation probability p for the Maxwell (thin continuous line), 
VHP (dashed line), and hard spheres models (thick continuous line). The vertical lines gives the value p = 0.893 . . . for which 
a divergence of the hard sphere transport coefficients ft* and fi* appears (while the shear viscosity exhibits regular behavior, 
see Fig. |3J . 
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FIG. 5: Transport coefficient n* as a function of the annihilation probability p (see Fig. |1] for more details). The Maxwell 
model is not represented since in this case fi* = 0. 

VII. STABILITY ANALYSIS OF THE NAVIER-STOKES HYDRODYNAMIC EQUATIONS 

A. Dispersion relations 



The hydro-dynamic equations. 135(1 cannot be solved analytically in general. However, their linear stability analysis 
allows one to answer the question of formation of spatial inhomogeneities. The present study establishes under which 
conditions the homogeneous state is stable. We consider here a small deviation from spatial homogeneity [see Eqs. i|18|) 
and l@2J] and the linearization of Eqs. l-iol) in the latter perturbation. The procedure used here follows the same 
route as for granular gases [^2 or PBA of hard spheres [28j . We define the deviations of the hydrodynamic fields 
from the homogeneous solution by 5y(r,t) — y(v,t) — yn{t), where y = {n, u, T}. Inserting this form in the Navier- 
Stokes-like equations yields differential equations with time-dependent coefficients. In order to obtain coefficients 
that do not depend on time, it is necessary to introduce the new dimensionless space and time scales defined by 
1 = v h (i)\/m/[fcBTfj(t)]r/2, r = f Q cLsvqh{s)/2, as well as the dimensionless Fourier fields Pkfj) = Srik(T) /nn (t) , 

w k (r) = y/m/[kBT H (T)]Su k (T), and k (r) = 5T±{t)/T h {t), where 5j/ k (r) = J Rd dl e-^ l 5y(l, r). Note that 1 is 
defined (up to a constant prefactor) in units of the mean free path for a homogeneous gas of density njj(t). The 
dimensionless time r(i) gives the accumulated number of collisions per particles up to time t. Since we will study 



both the Maxwell and VHP systems, we recall here the general results valid for non- vanishing decay rates 
and £u ■ Making use of the dimensionless variables, the linearized hydrodynamic equations read 



(0) 



AO) 

ST ' 



dr P ^ T 



d-1 



* 1 1 
-Tj k 



9t +P^t - 2 (d-l) 



d + 2 2 
k k 



dr + VKn 



w k|l + ik 



(1 - K>*)Pk(r) + (1 - <X)0 k (r) 
1 



— - i.^ ^ 
dr 



-r,*k 2 



d+2 
2{d-l)' 



L*k 2 



0. 


(58a) 


0. 


(58b) 


0. 


(58c) 


0. 


(58d) 
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where the transverse mode = Wk — appears to be decoupled from the other equations. The longitudinal 
velocity field is given by Wk» = (wk ■ ek)ek, and e k is the unit vector along the direction given by k. The transversal 
velocity field w k± consequently defines (d — 1) degenerated shear modes. Upon direct integration, we have 

w ki (t) = w k± (0) exp[s_L (p, k)r] , (59) 

where 

a±(p,k)= P $ ) *-~v'k 2 . (60) 

On the other hand, the longitudinal velocity field Wk, lies in the one dimensional vector space generated by k. Hence 
there are three hydrodynamic fields to be determined, namely the density p k , temperature #k, and longitudinal velocity 
field wt. = Wk|,<3k. The hydrodynamic matrix M of the corresponding linear system reads off straightforwardly from 
Eqs. Ij58b|) - (|58d|l . The corresponding eigenmodes are given by <p n (k) — exp[s„(p, k)r], n = 1, . . . , 3, where s„(p, k) are 
the eigenvalues of M. Each of these three fields is a linear combination of the eigenmodes; thus only the biggest real 
part of the eigenvalue s n (p, k) has to be taken into account to discuss the limit of marginal stability of the different 
modes. 

We define k± by the condition Re[s^(fc_L,p)] = 0, i.e., 

k x = [2i4 0) 7t7*] 1/a > (61) 

and A: 1 1 by max^ Re[s||(A;||,p)] = 0, k\\ < kx- Therefore if k > k± all rescaled modes are linearly stable. For k G 
only the rescaled shear mode is linearly unstable (the latter may however be non- linearly coupled to the other modes), 
and for k < k» all eigenvalues are positive which leads to instabilities. However, it should be kept in mind that 
the previous discussion involves rescaled modes only, and should be connected to the original r variable. Indeed, 
for any real system (for example a cubic box of volume L d ) the smallest wavenumber allowed for a perturbation is 
given by 2tt/L, which corresponds to the time-dependent dimensionless wavenumber fc m j n = 2ir/ (Ltwt C) where 
C = Att^ 1 ^ 2 /[{d + 2)T(d/2)]. Since the density n(t) is a decreasing function of time, fc m ; n increases monotonously 
and there exists a time t±_ such that k min (t) > fcj_ for t > tx- The lower cut-off k m - ln therefore eventually enters the 
region where the homogeneous solution is stable. For t = tx, the system is however not in a spatially homogeneous 
state, but it is nevertheless tempting to conclude that the perturbations will be damped for t > t±. Although this 
statement is not rigorously derived, we conclude here that an instability can only be a transient effect |3jj . 

The time tx can be estimated from the condition k m i n (t±) = kx- Making use of the hypothesis of small spatial 
inhomogcncities, we may replace the density n(t) appearing in the definition of k m i n (t) by the homogeneous density 
n#(i) given by Eq. (|47a|l . We obtain 

l/7n 

-lL (62) 

Is the transient instability alluded to easily observable in a simulation? A typical number of particles for molecular 
dynamics simulations is of the order of 10 5 , and n Q a 2 = 5 x 10~ 3 (which corresponds to a rather low total initial 
packing function 7m cT 2 /4 ~ 0.004). For p = 0.1 and d = 2 Eq. gives tx ~ 8.6 t Q . . .. Making use of Eq. I47a(l to 
approximate the density, one obtains n(t±) ss 0.61no. The density inhomogcncities therefore start to decrease after 
that the density decreased to only 0.61 times its initial value, which for p = 0.1 corresponds in average to only 4 
collisions per particle. For comparison purposes, inhomogeneities in granular gases are observed after a few hundred 
collisions per particle UH^]. In order to observe the previous (and presumably transient) instabilities one would 
need molecular dynamics simulations with very large systems. Another condition is to have a large enough p, which 
increases k±, see Fig. Equivalently, increasing p increases the divergence rate sj_ at fixed fc, see Eq. (|60|1 . For 
sufficiently small p (or small system sizes) Eq. (I62H does not have a positive solution because fc m ; n > k± already 
for t = 0. To sum up, the typical size of the inhomogeneities may grow as a function of time until t ~ tj_ but the 
subsequent evolution should drive the system back to a time dependent homogeneous regime. 

B. The linear second order contributions 

As already pointed out in Sec. IIVDI it is necessary for consistency to include the second order decay rates in the 
first order Navier-Stokes equations. Their study is useful to establish the relevance of the second order contributions 




(d + 2)T{d/2) 



kx(p) 
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in the gradients [that can be linear in V 2 n or V 2 T, or non-linear in (VT) 2 , VT ■ Vn, etc.]. Our analysis follows the 
method of Ref. 32], and we shall therefore not enter into much details. 

Since we shall perform a linear stability analysis, the only terms in the second order decay rates that will contribute 
are linear in the gradients. We therefore denote /j 2 ' the part of the second order distribution /( 2 ) that yields the linear 
contributions, and neglect any other term in . The solution has the form = M [n, T, Y)'S/ 2 T + N(n, T, V)V 2 n, 
where M and N are to be determined. Inserting the latter relation in the linear second order decay rates [that 



may be obtained from Eqs. I|21|l upon replacing by /j 2 '] they take the form £^ — s , , 



^\\7 2 n, where 



A = {n,Ui,T}. Again, it is useful to resort to a first order Sonine polynomial expansion for M and N, such that 



M(V) = c { rf> S 2 (c 2 )M(V) and N(V) = c\t' S 2 {c 2 )M{V), where ti£> and cit' are the coefficients to be determined. 

While it is a rather straightforward task to see from Eqs. lj2"Tl) that those second order contributions to the decay 
rates are equal to zero for the Maxwell gas, the calculations for the VHP model are more involved. In the latter 



case all decay rates are equal to zero except & 2 | = Xc, 



.-.(2) 



and £y 2 



Xc (2) 



coefficients c. 



r and Cn^ can be put in a dimensionless form <£p* = nfcs Tvq/kq and c^* = Cn'^ksvo/ Ko- Upon 



(2) 



where X = (T d - x <jmv T d{d + 2)/2. The 

„(2)^2, 



replacing the above expansions for and for the decay rates in the Boltzmann equation, the coefficients c^'" and 



Cn^* are solution of the equations 

2pd 0) 



(2)* 
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V 



(n)* 



where Y = Wtt^- 1 ^ 2 /[c/>r(d/2) y/2(d + 2)(d - 



2(d+2) 

d + 4 
2(d + 2) 
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d(d 



-K*Y, 



- p - d + 4 
4)], and 



.(2)* _ 



(0)* (2)* 



d(d + 2) 



2) 

H*Y, 



1 J Rd dVV 4 J[S 2 (c 2 )M(V)} A V2T(d/2) 



vo f ni dVV*S 2 (c 2 )M{V) 



47r( d - 1 )/ 2 



\p{d + 2)(d + 3) + (1 -p)4d] 



(63a) 
(63b) 

(64) 



The second order linear decay rates are given 

by^ 2 !* = c ( T } *d(d + 2)/[Y(d-l)(d + 4)] and^ 2 2* = c [ n ] *d(d + 2)/[Y(d- 

l)(d + 4)]. 

The dimensionless linearized hydrodynamic equations with nonzero second order decay rates are the same as 



Eqs. (|58(l except for Eq. (|58c() where one has to replace n* by k* — p£$p* an d M* by /i* 



P&1* 



Consequently, the 



ratios k* /pQp * an d M*/pCt2 ( see Fig.|HJ give a rough estimation of the relevance of the second order decay rates. 

As seen from Fig. |6l the second order decay rates [ and 2 are much smaller than the transport coefficients 
k* and /it*, respectively. They may therefore be neglected at least from a linear stability analysis point of view. This 
approximation is increasingly more accurate as the annihilation probability is decreased. Note that our conclusions 
for the Maxwell model of probabilistic ballistic annihilation are comparable to those for the granular gas [32|. 



C. Comparison between Maxwell, very hard particles and hard sphere results 

For the Maxwell model, the temperature decay rate vanishes. It follows from Eq. (|60(l that k± = and the 
transverse mode is stable, which is confirmed by Fig. [7| The Maxwell model appears to be linearly stable for all 
values of the annihilation probability p. On the other hand, within the VHP approach, the decay rate £t(0) ^ 0. The 
transverse mode may consequently be unstable for some wave- numbers k of the perturbation (see Fig.|S|), which by 
nonlinear coupling to the other modes may lead to density inhomogeneities. Other modes than the shear may also be 
linearly unstable, when rescaled wave numbers are such that k > fen . The thresholds k±_ and are shown in Fig. |5| 
for the 3 models. It appears again that the hard sphere quantity is bounded below by its Maxwell counterpart and 
above by VHP. Note that the linear stability analysis does not suffer from arbitrariness related to the free parameter 

The imaginary part of the eigenvalues embodies the information on the propagation of the perturbations. In Fig. EI 
we identify 3 different parallel modes for small enough k (k < 0.05). Given that the shear mode is always (d — 1) 
times degenerated and that there are d + 2 modes in total, none of the parallel modes are degenerated for low enough 
k. Increasing k up to the first bifurcation, the sound modes become degenerated and have a nonzero imaginary value. 
The non-propagating sound modes thus have bifurcated into a pair of propagating modes. Since the eigenvalue for the 
transverse velocity field is always real, we shall study here only the imaginary part of the other eigenvalues. We define 
k p such that for all k < k p all eigenvalues are real. It means that only perturbations with small enough wave numbers 
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FIG. 6: Ratio of the dimensionless transport coefficient k* (/i*) to the linear second order dimensionless transport coefficients 
£t i* f° r the VHP model. These ratios do not depend on the the collision frequency <f>- 
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FIG. 7: Real part of the eigenvalues in dimensionless units for the Maxwell model with p — 0.1 and d — 3. The dispersion 
relation obtained from Eq. 160H is represented by a dashed line (labeled sj_) whereas the three remaining relations are represented 
by continuous lines (labeled sn). The shear mode (s_l) and sound modes (which are on this figure such that s = when k — > 0) 
are degenerated twice. 
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FIG. 8: Real part of the eigenvalues in dimensionless units for the VHP model with p = 0.1 and d = 3. The dispersion relation 
obtained from Eq. 1601 is represented by a dashed line (labeled s±) whereas the three remaining relations are represented by 
continuous lines (labeled s«). The first two biggest parallel modes are sound modes. 




FIG. 9: Wavenumber k± and fey in dimensionless units as a function of the annihilation probability p for d — 3. S and 
VHP superscripts denote the hard spheres and very hard particles models, respectively. Within the Maxwell model, one has 
k± = fey = 0. 
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FIG. 10: Wave number k v in dimensionless units as a function of the annihilation probability p for d — 3. The Maxwell model is 
not represented since in this case k p — for all p. The main inset shows the imaginary part of the eigenvalues in dimensionless 
units for the VHP model for d = 3 and p = 0.1. The smaller inset shows the propagation gap k £ [0.1006 . . . , 0.1046 . . .] of the 
sound modes. 



A such that A -1 > k p / (2irna d ~ 1 ) are propagating. Fig. I1UI shows k p as a function of the annihilation probability p for 
the VHP, hard sphere, and Maxwell models. Once more, the VHP and Maxwell models appear as upper and lower 
bounds, respectively. From Fig. the Maxwell sound modes are degenerated for all k and therefore the sound modes 
of the Maxwell model are always propagating, i.e., k p = 0. In the VHP case, Fig. [H] shows a propagation gap for 
the sound modes, i.e., a k- window with k > k p , where the sound modes are not degenerated. This is confirmed by 
Fig. ITUI fsmaller inset). A propagation gap in the sound mode dispersion relation has been predicted on the basis of 
the revised Enskog theory for hard-sphere fluids [3^ . Such a gap has been observed in neutron scattering experiments 
of atomic liquids |4(j, of molten salts |4l|, or inelastic x-ray scattering of lipid bilayers 01 f° r example. 



VIII. CONCLUSION 



Making use of the Chapman-Enskog scheme, we have derived in this paper the hydrodynamic equations governing 
the coarse-grained number density, linear momentum and kinetic energy density fields for an assembly of particles 
undergoing annihilating collisions with probability p and an elastic scattering otherwise. In between collisions, the 
motion is ballistic. To this aim, the relevant "hard sphere" -like Boltzmann equation has been simplified first into its 
Maxwell, and second into its very hard particle (VHP) form. In both cases, the corresponding Navier-Stokes equations 
take the same form as in the initial hard sphere description and read 

d t n + V • (nu) = — pn £„, (65a) 

d t u+— V-P + (u- V)u = -pv T £ U7 (65b) 
mn 

dtT + (u • V)T + -t^(P : Vu + V • q) = -pTfr, (65c) 
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with 



= -^ 4 TP) M V^T' (66) 

£u = 0, (67) 
Ct - 0, (68) 



(69) 



e« = 


2d 




d + 4 




€u = 


-v T ^ 


K* 


= 


2 




d + 4 


l/ 0> 



for the Maxwell model, and 



T " f ^ 7i "7 2(d-l)(d + 4)' (?0) 

(71) 

in the VHP case [the transport coefficients k* or fx* are given by Eqs. (|50|l ]. 

Our analysis showed that the Maxwell and VHP simplifications, that are more amenable to analytic treatment, not 
only capture the essential features of hard sphere dynamics, but also provide lower and upper bounds for all comparable 
quantities. Some important differences should however be commented upon. A first difference is that Maxwell and 
VHP lead to regular transport coefficients for all values of the annihilation probability whereas a divergence occurs 
for annihilating hard sphere thermal conductivity K and Fourier coefficient /i. We concluded from this comparison 
that this divergence is presumably not physical and could result from the more stringent approximations put forward 
in the hard sphere computation. It turns out that the hard sphere case is such that the velocity distribution is 
non-Gaussian to zeroth order in spatial gradient, whereas it is Gaussian in Maxwell and VHP cases. This fact could 
be at the root of the divergence observed in the transport coefficients. We have also shown that the second order 
decay rates of the Maxwell model are equal to zero, while for the VHP case they may accurately be neglected. This 
analysis suggests that the second order decay rates of probabilistic ballistic annihilation of hard spheres may as well 
be accurately neglected, at least for small annihilation probabilities p. This approximation had been invoked in Ref. 
|28j , without any control on its validity. 

The second important difference between Maxwell, hard sphere and VHP dynamics is that within the Maxwell 
model, all Fourier modes are found to be linearly stable. This fact is intimately related to the non dissipative nature 
of the corresponding dynamics, an aspect which may be surprising at first: although particles are permanently removed 
from the system, the mean kinetic energy is conserved on average (£t = 0). This may be considered as a deficiency 
of the Maxwell (over)simplification. On the other hand, VHP dynamics is such that the collision frequency increases 
with the velocity of a given population of particles, which in turn implies that the kinetic energy decreases faster than 
the number of particles, hence £t > 0. This dissipation is at the root of possible instabilities in the coarse-grained 
fields. However, these instabilities manifest themselves for suitably rescaled fields, and we argued in section IVTT1 that 
they should presumably only translate into transient instabilities for the "real" fields. Indeed, due to the decrease of 
density n(t), an unstable Fourier mode has a wavenumber increasing like and eventually enters into a regime 
where damping should wash out the perturbation. This feature presumably provides at least a linear saturation 
mechanism for instabilities, different from usual non-linear saturation effects, that may also play a transient role here 
if the initial conditions are sufficiently unstable [in other words, if n(t±) <C n(t = 0)]. Our stability analysis was 
nevertheless restricted to perturbations around the time dependent homogeneous state, so that strictly speaking, the 
effects of transient instabilities that may drive the system into a strongly modulated state are unclear at the moment. 
This calls for a careful numerical (molecular dynamics) study of the coarse-grained fields, which is the purpose of 
future work. This would also allow to question the validity of the hydrodynamic description, in a regime where the 
wave number is not much smaller than the inverse mean-free-path l~ x cx ncr d_1 (in the previous Figures, k is expressed 
in units of t~ x , up to a prefactor of order one). 
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APPENDIX A: SUMMARY OF THE NOTATIONS 

We shall recall here some of the notations used through the paper, k and [i are the transport coefficients appearing 
in Fourier's linear heat conduction law Q24[l. and r\ is the shear viscosity appearing in the pressure tensor i|23|) . A 
quantity A that is made dimensionless is noted A*. The corresponding dimensionless transport coefficients are written 

»?* = -, (Ala) 



n* = — , (Alb) 
= (Ale) 



where 

2(d- 1)^''°' 



d(d + 2) fc B 

K = 7773 TT %, (A2) 



d + 2 r(d/2) 

are the thermal conductivity and shear viscosity coefficients for hard spheres, respectively. The dimensionless coeffi- 
cients u* v*, and v* are given by 



1 J Rd dV S l (V)JA l 1 J Rd dV Si(V)nAi 



v J Rd dVS t (V)A z 1 w J Rd dVS t (V)A t ' 



(A4a) 



P.. r Jir o r\r\K> ' (A4DJ 



. 1 / Rd dVA 3 (V)JC I3 1 J KJ rfVA 3 (V)^ 



p (0) g 



with 

= - — = , - - '' , rr-r na d ~ L \/-^, (A5) 
% d + 2 r(d/2) V m ' v ; 

and p^ - 1 = nksT is the zeroth order pressure. In Eqs. I|A4(1 . the operator J is given by 

Jg = P L a [fW,g] + (l-p)L c [f i0 \g], (A6) 

where 

L a [f i0 \g] = -Ja[f {a) ,g}-Ja[g,f { % (A7a) 
L c [f {0 \g} = ~Jc[f i0 \g}~J c [gJ { % (A7b) 

(A7c) 

g being an arbitrary function. The collision operator J c (annihilation operator J a ) is defined by Eq. © [Eq. (J2J)]. 
The linear operator is defined by Eq. (|D3jl . 

The velocity distribution function is denoted /(r, v;t). In the scaling regime 

'< r ' v;i > = H M (A8 » 

where c = V/vt- The time dependent [through T(t)] thermal velocity is 



«t = V — — . (A9) 

TO 
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We note the Maxwellian in the homogeneous cooling state by 

"(*) 
v d (t)n d / 

and the Maxwellian by 

M(c)=7T- d/2 cxp(-c 2 ). (All) 



Therefore, we obtain a similar relation to Eq. I|A8|) : M(V) — (n/v^)A4(c) 

The decay rate for the field A = {n, m, T} re< 
The corresponding dimensionless decay rate is 



The decay rate for the field A — {n, m, T} reads , where m denotes the order in the Chapman-Enskog expansion. 



»(to) 
Ma 



gm)* = U_ (A12) 



APPENDIX B: DECAY EXPONENTS OF THE HOMOGENEOUS COOLING STATE OF PBA WITH 

HARD SPHERE DYNAMICS 



We shall recall results obtained in [271 |28j . The density and temperature of the homogeneous cooling state for PBA 
of hard spheres are given by 

n H {t) = n fl+pl^ , (Bla) 

T H (t) = T U+P^J T > ( Blb ) 

where the decay exponents are 7^ = ^i°^(0)/t o , 77 = 6^(0)/io, and the relaxation time to = £n\o) + £^(0)/2, 

where £n\o) and £^(0) are the decay rates at time t — 0. The subscript H denotes a quantity evaluated in the 

homogeneous state. Making use of the explicit values of ^''(O) and £^(0) found in j2^| and since 7^ = 7j./2 one 
obtains 

q d + 2 ( 1 \ , 

7^ = —r- 1 ~ a 2-^ v to, (B2a) 



4 V 16, 
q ld + 2 ( 8d+ll\ 



to 1 



d + 2 ( 1 \ d + 2 ( 8d + 11 

1 - 02-7: + — TT-T 1 + 0,2 



4 V 16/ 16d \ 16 

where the kurtosis a 2 of the velocity distribution is [2 



vt>, (B2c) 



a 2 = 8 = ^ = . (B3) 

4d + 6-x/2 + i^872(d-l) V ; 

APPENDIX C: USEFUL RELATIONS FOR THE COEFFICIENTS v% AND v* 

The expressions (12811 and (|51|l may be calculated with the help of the following relations. Let X and Y be arbitrary 
functions, A4(V) = n/(v^TT d ^ 2 ) exp(— V" 1 jv\~) the Maxwellian in the scaling regime, then 

/ dv 1 Y(v 1 )L a [MX]=a d - 1 <p(x)v^- x [ dw l dv 2 < 2 / (0) (vi)M(v 2 )X(v 2 ) [F(vi) + Y(y 2 )} , (CI) 



and 



dv 1 Y(vi)L c [MX] = _ a d-i 4>{x)v l T x I dv 1 dv 2 vf 2 f(°\v 1 )M(v 2 )X(v 2 ) [ da(b - 1) [F( Vl ) + Y(v 2 )] , (C2) 
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where L a g = -Ja[/ (0) ;.g] - J a [g,f^] and L a g = -./ c [/ (0) ,g] - J c [g, / (0) ] for an arbitrary function g. Let a S M+, 
then 

ixlxl-e— 2 - ^ r[(d + »)/2] 

5(X|X| 6 _ a(d+n)/2 r(d/2) ' [ ^ 6) 
, I m -ax 2 *" d/2 d + nrlfd + n)/^ 

T1 / x|x| e ^ = swid — fW2r <5 - (C4) 

In the integrals below, the results when 6{a ■ g) is absent are obtained upon multiplying the value of f3 n by two. 
Finally for a = (oi, . . . , a d ), g € R d , |tr| = 1, |g| = 1, we have 

/" 9(9 ■ g)(? • g)"^- = -A^ 5 «- 2 (n 3 ^ + ff 2 %), (C5) 

/" 0(5 • g)(S ■ g) V, = /3„ +lff "- 1 5l , (C6) 

A = /^ e(S .B( S . ir = ^£|±|Z|. (C7) 

APPENDIX D: EXACT RELATIONS FOR THE TRANSPORT COEFFICIENTS OF THE MAXWELL 

MODEL 

Following the same route as in [2^ we may rewrite the right hand side of Eq. 1|19|) such that 

[^ 0) + J]/ (1) =AiVilnT+BiVilnn + CijV i u j +pClf^, (Dl) 

where 

Ai = ^A W/ (0)] _ (D2a) 
2 9Fj 1 J,/ J to 0^ ' v ; 

« . - (D2b) 
to aVi 

C " - asM/^-i^W/Wft., (D2c) 

and f2 is a linear operator defined by 

= /<°>tf - + (D3) 

The quantity g is either B,-, or Cy, and the functionals Qt\ (,u) j and £^ are obtained from Eqs. (|21(l upon 
replacing /W by 3. 

1. Pressure tensor 

Integrating the Boltzmann equation (|D1(1 over V with weight mViVj and taking into account the symmetry prop- 
erties of the coefficients (|D2ll one obtains 

d^pV(r,t)+p [ dVmViVjLalfWj^ + il-p) [ dVmViVjL B \fW, = / dV mV^C^V)^ k u h (D4) 
jR d Jn d J~R d 

where we have made use of the definition (|25() for the pressure tensor. The same definition further allows us to write 

/ dVmV i V j L a [f i0) ,f (1) }=e ) P^(r,t), (D5) 
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and using additionally Eq. 1|C2|) , and l|C5f) to (JC7I) : 

Finally 



where 



Aijfe; — SikSji + SjkSu — —SijSki- 



Insertion of Eqs. l|D5|) to (|D7)l in (|D4|l yields 



d + 2 



P^\r,t)=-p^A ljkl V 



kUl- 



(D6) 



(D7) 



(D8) 



(D9) 



The solution of Eq. is P^' = _r ? A y 

feiVfeU;. Functional dependence analysis shows that 77 oc T 1 / 2 , and since to 
zeroth order the temperature is conserved d t P^ = 0. Eq. (|D9|I thus gives 

' (D1Q) 



P^ + (1-P) 



2. Heat flux 



Integrating the Boltzmann equation (|D1(1 over V with weight mV 2 Vi/2 and taking into account the symmetry 
properties of the coefficients (|D2|) one obtains 



[ dv\mV 2 V l A k {Y)V k \nT + / dV imI/ 2 V^ fe (V)V fe InT, (Dll) 



where we have made use of the definition for the heat flux to first order. Moreover, one finds 

/ dvlmV 2 V i L a [f(°),fU] = ^ ) qV(r,t), 

J]R d 1 

and using additionally Eq. JU2J), and to (fU7|l: 



Finally 



/ dv\mV 2 V t A k {\)V k \nT = - 



d(d + 2) 

d + 2p^k B 



2 m 



ViT, 



/ dvim^ 2 ^B fc (V)V fc lnn = 0. 

JTR d 2 

Insertion of Eqs. HD12|) to HD15|I in (|D1 1|) gives 



(0) 



2(rf- 1) 
d(d + 2) 



<z, (1) M) = - 



2 m~ 



(DI2) 



(D13) 



(D14) 



(DI5) 



(DI6) 
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The solution of Eq. (|D16|I is gj 1 -' = — AV^T — /iV;n. Functional dependence analysis shows that A oc T 1 / 2 and 
fi oc T 3 / 2 n -1 , therefore dtq^ = p£n '/A/ i n - I n Q1 'der to satisfy Eq. IJD16I) it is therefore required that /i = and 

A * = i h — ( D1? ) 

PW^Tj + (1-P) 



Granular Gases, edited by T. Poschel and S. Luding, Lecture Notes in Physics vol 564 (Springer- Verlag, Berlin, 2001). 
J.W. Dufty, Recent. Res. Devel. Stat. Phys. 2 , 21-52 (2002) (e-print: cond-mat/0108444 ); J.W. Dufty, Adv. Compl. Syst. 4, 
397-406 (2001) (e-print: cond- mat/0109215 ). 

N.V. Brilliantov and T. Poschel, Kinetic Theory of Granular Gases (Oxford University Press, 2004). 
J.W. Dufty and J.J. Brey, e-print: cond-mat/0410133 

M.H. Ernst and E.M. Hendriks, Phys. Lett. 70A, 183 (1979); E.M. Hendriks and M.H. Ernst, Physica 120A, 545 (1983). 
M.H. Ernst, Phys. Rep. 78, 1 (1981). 

U. Marini Bettolo Marconi and A. Puglisi, Phys. Rev. E 66, 011301 (2002). 
A. Baskaran and J.W. Dufty, e-print: cond-mat/0410084 
E. Ben-Nairn and J. Machta, |cond-mat/0411743| 
A. Santos, Physica A 321, 442 (2003). 

J.W. Dufty, A. Baskaran, and L. Zogaib, Phys. Rev. E 69, 051301 (2004). 

V. Garzo and A. Astillero, e-print: cond-mat/0404386 

V. Garzo and J.M. Montanero, e-print: cond-mat/0411221 

M.J.E. Richardson, J. Stat. Phys. 89, 777 (1997). 

N.N. Bogolubov, in Studies in Statistical Mechanics, edited by J. De Boer and G.E. Uhlenbeck, Vol. I (North-Holland 
Publ. Co., Amsterdam, 1962), p. 5. 

E.G.D. Cohen, in Statistical Mechanics of Equilibrium and Non-Equilibrium, edited by J. Meixner (North-Holland Publ. 
Co., Amsterdam, 1962), p. 140. 

S. Chapman and T.G. Cowling, The Mathematical Theory of Non-Uniform Gases (Cambridge University Press, Cambridge, 
1970). 

Y. Elskens and H.L. Frisch, Phys. Rev. A 31, 3812 (1985). 

E. Ben-Naim, S. Redner, and F. Leyvraz, Phys. Rev. Lett. 70, 1890 (1993); E. Ben-Nairn, P.L. Krapivsky, F. Leyvraz, and 
S. Redner, J. Phys. Chem. 98, 7284 (1994). 
J. Piasecki, Phys. Rev. E 51, 5535 (1995). 

P-A. Rey, M. Droz, and J. Piasecki, Phys. Rev. E 57, 138 (1998). 

R.A. Blythe, M.R. Evans, and Y. Kafri, Phys. Rev. Lett. 85, 3750 (2000). 

P.L. Krapivsky and C. Sire, Phys. Rev. Lett. 86, 2494 (2001). 

J. Piasecki, E. Trizac, and M. Droz, Phys. Rev. E 66, 066111 (2002). 

E. Trizac, Phys. Rev. Lett. 88, 160601 (2002). 

F. Coppex, M. Droz, J. Piasecki, E. Trizac, and P. Wittwer, Phys. Rev. E 67, 021103 (2003). 
F. Coppex, M. Droz, and E. Trizac, Phys. Rev. E 69, 011303 (2004). 

F. Coppex, M. Droz, and E. Trizac, Phys. Rev. E 70, 061102 (2004). 

E. Trizac and P.L. Krapivsky, Phys. Rev. Lett. 91, 218302 (2003). 
A. Santos and J.J. Brey, Phys. Fluids 29, 1750 (1986). 

R. Soto, M. Mareschal, and D. Risso, Phys. Rev. Lett. 83, 5003 (1999). 

J.J. Brey, J.W. Dufty, C.S. Kim, and A. Santos, Phys. Rev. E 58, 4638 (1998). 

J. Ferziger and H. Kaper, Mathematical Theory of Transport Process in Gases, North-Holland, Amsterdam (1972). 
V. Garzo and A. Santos, Kinetic Theory of Gases in Shear Flow (Kluver Academic Publisher, 2003). 

F. Coppex, M. Droz, J. Piasecki, and E. Trizac, Physica A 329, 114 (2003). 
J.M. Montanero and A. Santos, Granular Matter 2, 53 (2000). 

Transient instabilities were also predicted for viscoelastic granular gases with velocity-dependent restitution coefficient, see 
N. Brilliantov, C. Saluena, T. Schwager, and T. Poschel, Phys. Rev. Lett. 93, 134301 (2004). 
I. Goldhirsch and G. Zanetti, Phys. Rev. Lett. 70, 1619 (1993). 

M.J. Zuilhof, E.G.D. Cohen, and I.M. de Schepper, Phys. Lett. A 103, 120 (1984); E.G.D. Cohen, I.M. de Schepper, and 
M.J. Zuilhof, Physica B 127, 282 (1984); E.G.D. Cohen, Physica A 194, 229 (1993). 

A. A. van Well, P. Verkerk, L.A. de Graaf, J.-B. Suck, and J.R.D. Copley, Phys. Rev. A 31, 3391 (1985); A. A. van Well 

and L.A. de Graaf, Phys. Rev. A 32, 2396 (1985). 

R.L. McGreevy and E.W.J. Mitchell, J. Phys. C 18, 1163 (1985). 

T.M. Weiss, O.-J. Chen, E.E. Alp, S.-H. Chen, and H.W. Huang, Biophys. J. 84, 3767 (2003). 



